In [241]:
import numpy
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import *
Natural convection in cavity:
In this part we are going to simulate natural convection problem in a cavity. The cavity consists of four sides. The right wall is hot and opposite wall is cold. Top and bottom edges are considered insulated or adiabatic condition. Let’s have a look at the equations.
$$\frac{\partial(\rho \overrightarrow{u})}{\partial t} + \overrightarrow{\nabla}\cdot[\rho\overline{\overline{u\otimes u}}] = -\overrightarrow{\nabla p} + \overrightarrow{\nabla}\cdot\overline{\overline{\tau}} + \overrightarrow{F}$$$$\frac{\partial(\rho c_pT )}{\partial t} + \overrightarrow{\nabla}\cdot(\rho c_pT\overrightarrow{u}) = \nabla\cdot(\kappa(\nabla{T}))$$What are the added source terms (Fx and Fy ) on the right hand side of Navier stokes equations ? Thanks to the heat transfer course we know that in natural convection phenomena the gas is heated and subsequently the density of the gas is going to be changed. Therefore based on ideal gas hot gas has lower weight compared to cold gas. Hence the hot gas moves up and the cold gas moves in opposite direction (downward). This is the mechanism occurred near your room‘s radiator. Major changes are seen in these equations if we compare them with our previous examples. The first one is that velocity of the gas depends on temperature distribution and on the other side velocity of gas affects the temperature distribution. Therefore these equations are coupled. We need to solve them together and update the values in each time step.
As discussed earlier, for each problem we need to apply “chapman-enskog” relation to find that what the relation between relaxation time and diffusion parameter is. Let’s at first convert these equations to non-dimensional forms. You can review your heat transfer book to find how you can use characteristic and scales to make the equations non-dimensional.
$$\frac{\partial\overrightarrow{u}}{\partial t} + \overrightarrow{\nabla}\cdot[\overline{\overline{u\otimes u}}] = -\overrightarrow{\nabla p} +\nabla^2 \overline{u}+Ra\theta $$$$\frac{\partial \theta}{\partial t} + \overrightarrow{u}\cdot\nabla{\theta} = \frac{1}{Re\cdot Pr} (\nabla^2 \theta) $$Where:
If you review the lid driven code, we mentioned that it is necessary to make equal Reynolds number based on LBM units. Why, because you cannot use directly the values such as kinematic viscosity in LBM. Why? Because you have limitation in using Mach number and it will be instability problem if you use unreasonable value for your relaxation time or velocity speed characterization. On the other side, consider the cavity length as one. Can you use one node in your LBM code? Definitely your answer is not. You instead use for example 100 nodes to model the length. You are changing the number of nodes however you should be confined to the length one. Hence one relation should be made to convert the 100 nodes to the length of one. That is why we need to make equal dimensionless numbers. In FDM for having Reynolds number 100 you take the length of the cavity as one and the viscosity as 0.01 and the characteristic velocity speeds one. But here you are not allowed to choose one for your length instead we try to have the same dimensionless number which affect the physics of our problems. Using “chapman-enskog” specifies the Fx and Fy as follow:
$$F = 3\omega(k) g_x \beta \theta e_x +3\omega(k) g_y \beta \theta e_y $$That is why we need to have dimensionless equations. In this problem two dimensionless numbers affect the fluid behavior (Rayleigh and Prandtl). Increasing the Ra number leads to the natural convection and decreasing leads to conduction condition. Prandtl is also the criteria for measuring the hydrodynamic boundary layer to thermal boundary layer thickness.
Changes:
In this example we need to have two distribution functions, one for modelling fluid and the second for heat transfer (f and g). The streaming section does not change and it is similar to what we used in liddriven problem. We do not have source term in heat equation (g) but the source term should be added for fluid equation (f) . If you want to consider inclined cavity you need to use the complete formula for (F) however in this example as can be seen the force is in y direction and only the second part should be added. Each distribution function need each relaxation time based on diffusion parameter:
$$\omega_m = \frac{1}{3\nu+0.5} ~~~ (\nu=\frac{\mu}{\rho})$$$$\omega_s = \frac{1}{3\alpha+0.5}~~~(\alpha=\frac{k}{\rho c_p})$$
In [242]:
nx=50 # the number of nodes in x direction lattice direction
ny=50 # the number of nodes in y direction lattice direction
pr=0.71
visco=0.02 #*ny
dt=1.0 # temperature difference between left and right wall
alpha=visco/pr # heat diffusion coefficient
Tw=1.0 # left wall temperature
rhoo=1.0
D=3.0 # the dimension of the problem
omega=1./(0.5+(visco*D)) #Chapman-Enskog dt=1 and dx=1 relaxation flow
omegat=1./(0.5+(alpha*D)) # relaxation temperature
mstep=80000 # the number of time step
k=9 # k=0,1,2,3,4,5,6,8,9
In [243]:
w=numpy.ones(k) # witghting factor
cx=numpy.ones(k)
cy=numpy.ones(k)
rho=numpy.ones((ny+1,nx+1) ) # density matrix
u=numpy.ones((ny+1,nx+1) )
v=numpy.ones((ny+1,nx+1) )
f= numpy.ones((k, ny+1,nx+1)) # distribution function
feq= numpy.ones((k, ny+1,nx+1))
g= numpy.ones((k, ny+1,nx+1)) # distribution function for temperature
geq= numpy.ones((k, ny+1,nx+1))
T=numpy.ones((ny+1,nx+1) )
t1=numpy.ones((ny+1,nx+1) )
t2=numpy.ones((ny+1,nx+1) )
force=numpy.ones((ny+1,nx+1) )
In [244]:
strf=numpy.ones((ny+1,nx+1) ) # streamfunction
In [245]:
w[0]=4./9 #4./9
w[1:5]=1./9.
w[5:9]=1./36.
cx[0]=0.0
cx[1]=1.0
cx[2]=0.0
cx[3]=-1.0
cx[4]=0.0
cx[5]=1.0
cx[6]=-1.0
cx[7]=-1.0
cx[8]=1.0
cy[0]=0.0
cy[1]=0.0
cy[2]=1.0
cy[3]=0.0
cy[4]=-1.0
cy[5]=1.0
cy[6]=1.0
cy[7]=-1.0
cy[8]=-1.0
In [246]:
##================== Initial value
rho[0:ny+1,0:nx+1]=rhoo
v[0:ny+1,0:nx+1]=0.0
u[0:ny+1,0:nx+1]=0.0
T[0:ny+1,0:nx+1]=0.0
T[0:ny+1,0]=Tw
for i in range(nx+1):
for j in range(ny+1):
for l in range (k): #k=0,1,2,3,4
f[l,j,i]=w[l]*rho[j,i]
g[l,j,i]=w[l]*T[j,i]
Tref=0.5
In [247]:
## Main loop : comprised two parts :collision and streaming
##=====================
def function (Ra) :
gb=(Ra*visco*alpha)/(ny**3)
for n in range(mstep) :
if (n%10000==0):
print n
# collision process fluid flow
t1[0:ny+1,0:nx+1]=u[0:ny+1,0:nx+1]*u[0:ny+1,0:nx+1]+v[0:ny+1,0:nx+1]*v[0:ny+1,0:nx+1]
for l in range (k):
t2[0:ny+1,0:nx+1]=cx[l]*u[0:ny+1,0:nx+1]+cy[l]*v[0:ny+1,0:nx+1]
force[0:ny+1,0:nx+1]=(3.*w[l]*gb*(T[0:ny+1,0:nx+1]-Tref)*cy[l]*rho[0:ny+1,0:nx+1]*cos(alpha))+\
(3.*w[l]*gb*(T[0:ny+1,0:nx+1]-Tref)*cy[l]*rho[0:ny+1,0:nx+1]*sin(alpha))
feq[l,0:ny+1,0:nx+1]=w[l]*rho[0:ny+1,0:nx+1]*( 1.+ 3.*t2[0:ny+1,0:nx+1]+4.5*t2[0:ny+1,0:nx+1]**2-1.5*t1[0:ny+1,0:nx+1] )
force[0,0:nx+1]=0.
force[ny,0:nx+1]=0.
force[0:ny+1,0]=0.
force[0:ny+1,nx]=0.
f[l,0:ny+1,0:nx+1]=(1.-omega)*f[l,0:ny+1,0:nx+1]+omega*feq[l,0:ny+1,0:nx+1]+force[0:ny+1,0:nx+1]
f[1,0:ny+1,nx:0:-1]=f[1,0:ny+1,nx-1::-1]
f[3,0:ny+1,0:nx]=f[3,0:ny+1,1:nx+1]
f[2,ny:0:-1,0:nx+1]=f[2,ny-1::-1,0:nx+1]
f[5,ny:0:-1,nx:0:-1]=f[5,ny-1::-1,nx-1::-1]
f[6,ny:0:-1,0:nx]=f[6,ny-1::-1,1:nx+1]
f[4,0:ny,0:nx+1]=f[4,1:ny+1,0:nx+1]
f[8,0:ny,nx:0:-1]=f[8,1:ny+1,nx-1::-1]
f[7,0:ny,0:nx]=f[7,1:ny+1,1:nx+1]
## =============================
#left Boundary- Bounce back
f[1,0:ny+1,0]=f[3,0:ny+1,0]
f[5,0:ny+1,0]=f[7,0:ny+1,0]
f[8,0:ny+1,0]=f[6,0:ny+1,0]
#right Boundary-bounce back
f[3,0:ny+1,nx]=f[1,0:ny+1,nx]
f[7,0:ny+1,nx]=f[5,0:ny+1,nx]
f[6,0:ny+1,nx]=f[8,0:ny+1,nx]
# top Boundary -bounce back
f[4,ny,0:nx+1]=f[2,ny,0:nx+1]
f[7,ny,0:nx+1]=f[5,ny,0:nx+1]
f[8,ny,0:nx+1]=f[6,ny,0:nx+1]
# bottom Boundary -bunce back
f[2,0,0:nx+1]=f[4,0,0:nx+1]
f[5,0,0:nx+1]=f[7,0,0:nx+1]
f[6,0,0:nx+1]=f[8,0,0:nx+1]
for i in range(nx+1):
for j in range(ny+1):
sumr=0.0
for l in range (k):
sumr=sumr+f[l,j,i]
rho[j,i]=sumr
for i in range(1,nx):
for j in range(1,ny):
usum=0.0
vsum=0.0
for l in range (k):
usum=usum+f[l,j,i]*cx[l]
vsum=vsum+f[l,j,i]*cy[l]
u[j,i]=usum/rho[j,i]
v[j,i]=vsum/rho[j,i]
u[:,0]=0.
u[:,nx]=0.
u[0,:]=0.
u[ny,:]=0.
v[:,0]=0.
v[:,nx]=0.
v[0,:]=0.
v[ny,:]=0.
# collision process Temperature
for l in range (k):
t2[0:ny+1,0:nx+1]=cx[l]*u[0:ny+1,0:nx+1]+cy[l]*v[0:ny+1,0:nx+1]
geq[l,0:ny+1,0:nx+1]=w[l]*T[0:ny+1,0:nx+1]*( 1.+ 3.*t2[0:ny+1,0:nx+1] )
g[l,0:ny+1,0:nx+1]=(1.-omegat)*g[l,0:ny+1,0:nx+1]+omegat*geq[l,0:ny+1,0:nx+1]
#streaming process
# ==========================
# streaming Temperature
g[1,0:ny+1,nx:0:-1]=g[1,0:ny+1,nx-1::-1]
g[3,0:ny+1,0:nx]=g[3,0:ny+1,1:nx+1]
g[2,ny:0:-1,0:nx+1]=g[2,ny-1::-1,0:nx+1]
g[5,ny:0:-1,nx:0:-1]=g[5,ny-1::-1,nx-1::-1]
g[6,ny:0:-1,0:nx]=g[6,ny-1::-1,1:nx+1]
g[4,0:ny,0:nx+1]=g[4,1:ny+1,0:nx+1]
g[8,0:ny,nx:0:-1]=g[8,1:ny+1,nx-1::-1]
g[7,0:ny,0:nx]=g[7,1:ny+1,1:nx+1]
## Boundary condition
##============================ Temp boundary
#
#left Boundary
g[1,1:ny,0]=( Tw*(w[1]+w[3]) )-g[3,1:ny,0]
g[5,1:ny,0]=( Tw*(w[5]+w[7]) )-g[7,1:ny,0]
g[8,1:ny,0]=( Tw*(w[8]+w[6]) )-g[6,1:ny,0]
#right Boundary
g[3,1:ny,nx]=-g[1,1:ny,nx]
g[7,1:ny,nx]=-g[5,1:ny,nx]
g[6,1:ny,nx]=-g[8,1:ny,nx]
# top Boundary
for l in range(k):
g[l,ny,0:nx+1]=g[l,ny-1,0:nx+1]
# bottom Boundary
for l in range(k):
g[l,0,0:nx+1]=g[l,1,0:nx+1]
#
#
##================================
#
for i in range(1,nx):
for j in range(1,ny):
sumt=0.0
for l in range (k):
sumt=sumt+g[l,j,i]
T[j,i]=sumt
T[0:ny+1,0]=Tw
T[0:ny+1,nx]=0.
T[0,1:nx]=T[1,1:nx]
T[ny,1:nx]=T[ny-1,1:nx]
return
In [248]:
def plot (Ra):
strf[0,0]=0.0
for i in range(nx+1):
rhoav=0.5*(rho[0,i-1]+rho[0,i])
if(i>0) :
strf[0,i]=strf[0,i-1]-rhoav*0.5*(v[0,i-1]+v[0,i] )
for j in range(1,ny+1):
rhom=0.5*(rho[j,i]+rho[j-1,i])
strf[j,i]=strf[j-1,i]+rhom*0.5*(u[j,i]+u[j-1,i])
x=numpy.linspace(0,nx,nx+1)
y=numpy.linspace(0,ny,ny+1)
print '************* '
print 'Data for Ra=',Ra
print '************* '
print ' '
plt.figure(figsize=(30,5))
ax = plt.subplot(151)
plt.contourf(x,y,T,30)
title('Temperature Contour',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
ax = plt.subplot(152)
plt.contour(x,y,strf,30)
title('Streamline',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
ax = plt.subplot(153)
plt.contour(x,y,T,30)
title('Isotherm',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
ax = plt.subplot(154)
plt.contourf(x,y,u,30)
title('U contour',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
ax = plt.subplot(155)
plt.contourf(x,y,v,30)
title('Vcontour',fontsize=20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()
return
In [249]:
function (1e3)
plot(1e3)
In [250]:
function (1e4)
plot(1e4)
In [251]:
function (1e5)
plot(1e5)